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I.  INTRODUCTION 


An  ongoing  effort  is  being  conducted  at  the  US  Army  Armament  Research  and  Develop¬ 
ment  Command,  Ballistic  Research  Laboratory  (USA  ARRADCOM/BRL)  to  improve  the 
methodologies  that  are  used  to  calculate  the  vulnerability  of  military  vehicles  to  ballistic  threats. 
The  Internal  Point  Burst  model1  whose  distinguishing  feature  is  its  ability  to  separate  the  debili¬ 
tating  effects  of  spall  from  those  of  the  main  penetrator  constitutes  the  most  advanced  metho¬ 
dology  being  developed  at  the  BRL  and  elsewhere.  In  such  a  methodology,  the  vulnerability  of 
a  target  vehicle  is  calculated  as  the  expected  value  of  vehicle  incapacitation  due  to  the  direct 
impacts  of  primary  peneuators  and  any  associated  metal  debris  fragments. 

In  a  simplified  version  of  a  more  detailed  stochastic  representation  of  the  point-burst  vul¬ 
nerability  process2,3,  the  probability  equations  defining  the  expected  value  of  incapacitation  per 
projectile  (either  primary  penetrator  or  secondary  metal-debris  fragment)  for  a  vehicle  com¬ 
ponent  exposed  to  a  ballistic  threat  are  formulated  as  definite  integrals  whose  basic  form  is 
given  by 

A  = 

.wi.bi  lr0,w0,b0,g)0(«,i,w,,b,)</r,r/w1</bIrfr0</w0</b0t  (1A) 

go  oe  oo  oo  oo  oo 

J//5'(r0,w0,b())dr0</w0rfb0=*l.  (IB) 

oo  oo  oo 

In  this  expression,  the  boldface  lenars  (r^tvo,^)  and  (r),w,,b])  are  the  values  associated  with 
the  random  variables  (Ro,Wc,B0)  and  (Ri,Wt,B| ),  respectively,  where  R  quantifies  the  location 
and  W  quantifies  the  direction  of  motion  of  a  projectile.  Correspondingly,  the  set  of  variables  B 
is  used  to  quantify  (characterize)  some,  but  not  all,  of  the  remaining  significant  features  of  a 
projectile.  The  subscript  0  is  used  to  identify  the  variables  associated  with  a  penetrator  at  its 
origin  and  the  subscript  1  is  used  to  identify  the  variables  associated  with  a  penetrator  at  impact 
with  a  critical  component  (Figure  1)  in  the  vehicle.  The  quantities  d r,  dm,  and  rfb  are  the 
infinitesimal  hypervolumes  associated  with  r,  w,  and  b,  respectively. 

The  quantity  5(r0,w0>b0)  is  a  continuous  function  used  to  give  the  probability  density  that 
the  random  variables  associated  with  some  projectile  at  its  source  will  have  the  values 
(ro,w0,b0).  The  function/(ri,w1,bl  I  r^w^bog)  is  the  unnormalized  conditional  probability  den¬ 
sity  that  a  projectile  created  as  (r^w^bo)  will  perforate  a  barrier  characterized  by  the  set  of  ran¬ 
dom  variables  G  having  a  set  of  values  g  and  impact  the  critical  component  as  (rj.wj.bi).  This 
latter  function  differs  from  the  conventional  definition  of  a  probability  density  function  (PDF) 
in  that  it  does  not  normalize  to  unity,  but  has  a  normalization  which  for  each  set  of  values 
(r^w^b^g)  lies  on  the  range  from  0  to  1,  that  is 

0  ^  lr0lwo,b0,g)</ri</wi<fb1  <  I.  (1C) 

OO  OO  OO 

This  quantity,  herein  identified  as  the  perforation  PDF,  is  continuous  and  finite  in  bj  for  those 

1  J.R.  Rapp,  and  F  T.  Drown,  "An  Assessment  of  Existing  Point-Burst  Models  of  Armored  Vehicle  Vulnerability,”  Ballis¬ 
tic  Research  Laboratory  Memorandum  Report  No.  02963,  October  1979,  (UNCLASSIFIED).  (AD  #B043965L) 

2W.a  Beverly,  *A  Detailed  Stochastic  Ballistic  Vulnerability  Model  for  Armored  Military  Vehicles,”  Journal  of  Ballis¬ 
tics,  Vol.  ru.  No.  3,  1979. 

3w.a  Beverly,  ”A  Stochastic  Representation  of  the  Vulnerability  of  a  Critical  Component  m  a  MBHary  Vehide  to  Metal 
Fragments,”  Submitted  to  the  Journal  of  Ballistics  for  Publication. 
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Figure  l.  The  Mathematical  Notation  in  a  Simple  Spall  Example 


cases  where  a  birth  projectile  must  perforate  a  barrier  before  impacting  a  critical  component  in 
the  vehicle,  and  is  the  product  of  Dirac  delta  functionals4  in  bi  for  those  cases  where  the  pro¬ 
jectile  can  travel  unimpeded  to  impact  the  critical  component.  In  either  case,  the  perforation 
PDF  is  Dirac  delta  fUncuonu.j  in  r,  for  projectiles  impacting  the  critical  component  and  is 
identically  aero  over  its  total  range  for  projectiles  which  either  fail  to  perforate  the  barrier  or 
impact  the  critical  component.  The  function  Q(r,, is  any  well-behaved  function  whose 
values  on  the  impacted  surface  of  a  critical  component  give  the  avenge  incapacitation  caused  by 
all  penetrators  whose  "used*  random  variables  have  the  values  (ri,wltbi ).  In  these  integrals, 
S(ro.wOiho) IS  assumed  to  vanish  at  large  distances  from  its  maximum  value  so  that  an  integra¬ 
tion  over  all  values  of  all  variables  (identified  by  placing  the  infinity  symbol  at  the  bottom  of 
the  integral  sign)  will  yield  a  finite  expected  value  A  for  the  incapacitation. 

>  The  dimensionality  of  these  integrals  and  the  complexity  of  the  geometrical  and  composi¬ 
tional  features  of  the  target  vehicles  from  which  these  integrals  are  derived  could  make  their 
evaluation  by  deterministic  methods  prohibitively  expensive.  An  alternate  Monte  Carlo  method 
has  been  outlined  by  Beverly*  which  could  greatly  increase  the  efficiency  of  vulnerability  calcu¬ 
lations.  The  objective  of  this  study  is  to  analyze  the  Monte  Carlo  procedures  used  in  Reference 
S  and  to  illustrate  their  use  by  evaluating  simple  definite  integrals. 

In  the  next  section,  we  will  initially  analyze  the  Monte  Carlo  evaluation  of  simple  one¬ 
dimensional  integrals.  We  will  then  extend  the  procedures  developed  for  the  one-dimensional 
case  to  multiple  integrals.  We  will  also  analyze  integrals  having  the  form  illustrated  in  equation 
1A.  Then,  in  Section  ITT,  simple  integrals  will  be  constructed  which  have  the  form  of  the 
integrals  discussed  in  Section  IT.  These  integrals  will  be  solved  using  the  outlined  Monte  Carlo 
procedures  and  the  results  will  be  compared  with  the  known  closed  form  solutions. 


IT.  THE  EVALUATION  OF  DEFINITE  INTEGRALS  BY  USING 
THE  MONTE  CARLO  METHOD 


a.  A^.uBpip.QQ^Qimgnawoal.lflBBgraI 

A  simple  one-dimensional  integral  having  the  form 

A-jW)G(x)<fic,  (2) 

*1 

where  H(x)  and  G Or)  are  well-behaved  functions  on  die  interval  from  xj  to  xj  ,  is  generally 
regarded  as  the  area  under  the  curve  H(x)G(jc )  from  x,  to  x2.  However,  the  integral  can  be 
viewed  from  a  different  perspective  if  the  integrand  is  rearranged  to  obtain 

X-cjfmx)G(x)  dx-efHUnfbc)dx,  (3A) 

X,  U  — 


where 


C 


— fc(x)dx. 


(3B) 


4C  Eisenhan.  and  M.  Zelen,  "Elements  of  Probability,"  Handbook  of  Physics,  EU.  Condon,  and  H.  Odishaw,  Editors. 
McGraw-Hill  Book  Company,  Inc.,  New  York,  I9S8. 

5w.&  Beverly,  "The  Application  of  the  Monte  Carlo  Method  to  the  Solution  of  the  Internal  Point  Burst  Vehicle  Ballis¬ 
tic  VjJnerability  Model,"  Ballistic  Research  Laboratory  Technical  Report  in  Preparation. 
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and 

(3C) 

when 

x,  <x  <x2. 

(3D) 

and 

/(x)-0. 

(3E) 

when 

X  <X|,  x  >x2. 

(3H) 

The  integral  can  now  be  regarded  as  C  times  the  expected  value  of  Hix)  where  the  probability 

density  of  values  for  x  is  given  by  the  PDF/(x)  (Reference  4). 


Viewing  the  integral  from  the  second  perspective  and  applying  the  strong  law  of  large 
numbers  (Reference  4),  the  integral  can  be  estimated  as 

X— 7  (4a) 

J  j~\ 


where  the  xy  are  a  series  of  mutually  independent  values  for  x  whose  common  PDF  is  fix ). 
This  mean  value  estimate  is  identified  by  attaching  a  bar_to  A.  A  measure  of  the  confidence 
level  of  an  estimate  is  taken  to  be  its  standard  deviation  8A,  that  is6 


8A- 


icWXyWl* 

J-l 


/U-l) 


i 

I 


(4B) 


According  to  the  central  limit  theorem  (Reference  4),  this  measure  of  confidence  can  be  inter¬ 
preted  as  predicting  that  approximately  68  percent  of  a  large  number  of  similar  estimates  of  X 
will  fall  within  ±8A  of  X. 


A  step-by-step  outline  of  a  Monte  Carlo  evaluation  of  X  is  given  below  and  illustrated  in 
Figure  2. 

1.  Pick  a  value  x}  by  sampling  the  PDF  fix).  A  variety  of  methods  have  been  developed 
for  conducting  such  sampling7’*.  An  efficient  method  for  the  case  where  fix)  is  given  in  histo- 
gramic  form  is  to  pick  each  xy  by  solving  the  integral  equation9 

Jf  (x)dx  -*[iW(0, 1  ]j .  (5) 

The  quantities  £M0,l]y  are  a  set  of  independent  random  numbers  where  each  random  number 
is  picked  with  equal  probability  on  the  range  of  0  to  1.  The  reader  should  note  that  the 


*Y.  Beers,  "Introduction  to  the  Theory  of  Error,*  Addison- Weeley  Publish  mi  Company,  Inc.,  Hauling,  Mess.,  1962. 

7w.  Ojber,  R.  Nagel,  R.  Goldstein,  P.S.  Mittelman,  and  MJt  Katos,  *A  Geometric  Description  Technique  for  Com¬ 
puter  Analysis  of  Both  the  Nudeer  and  Convent  tonal  Vulnerability  of  Armored  Military  Vehicles,"  Mathematical  Appli¬ 
cations  Group,  Inc.  Report  No.  MAGI -6701 ,  August  1967. 

*E.D.  Cashwell,  GJ.  Everett,  and  O.W.  Rechard,  *A  Practical  Manual  on  the  Monte  Carlo  Method  for  Random  Walk 
Problems,’  Los  Alamos  Scientific  Laboratory  Report  No.  LA-2120,  December  1957. 

*Y.A  Schneider,  "The  Monte  Carlo  Method,*  Pergamon  Press,  Long  bland  City.  New  York.  1966. 
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Figure  2.  The  Mo  me  Carlo  Estimation  of  a  Simple  One-Dime  nsiomi  Integral 
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sampling  as  described  by  equation  S  generally  requites  that  the  PDP  be  normalized  to  unity.  We 
will  assume  that  this  method  is  used  in  the  following  outlines  of  Monte  Carlo  procedures. 

2.  Calculate  the  event  score  Xj  as  C  Hixj). 

3.  Accumulate  the  score  in  a  bin  reserved  for  this  operation. 

4.  Calculate  the  square  of  the  preceding  score. 

5.  Accumulate  the  squared  score  in  another  bin  reserved  for  this  operatioa 

6.  Reiterate  steps  1-5  for  a  total  of  J  everus. 

7.  Calculate  A  using  equation  4A. 

8.  Calculate  the  standard  deviation  8A  using  equation  4R. 

9.  Determine  if  the  confidence  level  of  A  as  determined  by  8A  is  adequate.  Conduct  more 
events  and  merge  their  results  with  those  already  obtained  if  8A  is  too  large. 

Step  9  completes  an  outline  of  a  Monte  Carlo  evaluation  of  a  one-dimensional  integral. 
The  procedures  outlined  above  will  be  extended  in  the  following  Subsection  to  the  evaluation 
of  multiple  integrals. 


B.  Multiple  Integrals 


A  multiple  integral  of  dimensionality  /  which  has  the  form 

A-/*  “JffU i, ...  ,Xi . X/)* 

A,  A, 

GCxj, .  . .  ,Xi, . . .  , Xj)dx\  •••<*}•••  dxj, 


(6A) 


can  be  evaluated  by  reiterating  the  procedures  outlined  in  Section  1TA  for  evaluating  one¬ 
dimensional  integrals.  The  integrations  are  assumed  to  be  started  at  the  innermost  integral  and 
conducted  toward  the  outermost  integral  where  1  is  the  running  index  over  the  variables.  The 
A,  and  B,  define  the  boundaries  of  the  integration  region  accordingly,  that  is 


A  i  —  constant, 
42=*  A2(X)), 

•  • 


A,  ^A1(xl. . 


B  j  =• constant ; 
fij”Bj(xt>, 


•  • 

Bf 


•  • 


A i  'mA/(x\,  ■ 


» xT _j  I,  Bj  —B/(X) i ...  i  Xj i ...  i xf  —  j )■ 


(6B) 


We  will  first  rearrange  the  integrand  in  equation  6A  to  obtain  : 
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X-C,  /•••/•••  Jffixu - X, . *f)* 

fix  . x, . Xf)d t\  •••<*;••  tbtu 


where 


Ci-/*  •  *  S'  **  /<3(x„  . .  •  ,x>. .  .  .  ,x,)<*r,  ••  •  dx,  -  dxr, 

A\  A,  A, 

GUf  .  ■ .  ,Xi,. . .  ,xr ) 


and 


/(*!,  .  .  .  ,Xj,  .  . .  ,X,)* 


(7A> 

(7B) 


(70 


when  the  x;-  have  values  located  within  the  region  bounded  by  the  A,  and  B,  ,  and 

fix i, .  .  .  ,Xj . X;)— 0,  (7D) 

when  the  x;  have  values  which  lie  outside  the  region  bounded  by  the  A,  and  B,.  Then,  using 
the  same  perspective  as  that  used  in  Subsection  IT  A,  the  multiple  integral  of  equation  7B  can  be 
interpreted  as  the  expected  value  of  Hixt, . . .  ,xt, . .  .  ,xr)  where  the  variables  x;  have  values 
predicted  by  the  joint  PDF  f(xu  .  .  .  .  xi . X/). 


However,  in  a  Monte  Carlo  estimation  of  X,  a  sample  value  (x)  )j  for  each  x  during  trial  j 
has  to  be  picked  individually  and  in  its  turn  The  probability  density  of  values  for  x,  is  the  mar¬ 
ginal  PDF  fiix\)  which  is  obtained  by  the  integration 

a,  B, 

f]bc i)— J  •*•/•••  ffixi. . .  ■  ,Xj, .  .  .  ,xj)dxi  •  •  •  dxi  •  ■  •  dxj.  (8a) 

A1  A,  A, 


Now,  similar  to  the  procedure  used  in  Section  TIA,  a  sample  value  (x,)y  can  be  derived  from 
the  integral  equation 


•*i 


Continuing,  a  marginal  PDF/2[(xi)/,x2]  is  then  constructed  for  the  variable  x2  as 
/2[(x,)y,x2]- 
^  Bl  Bi  Br 

7=-J  * '  *  J  •  •*  •  .x, . x,]<fcj, . . .  ,dxt,. . .  ,dxu 


(80 


where 

C2-J  *  *  ’  /  •  *  *  )Aix ,);,x2 . X),  .  .  .  ,xi\dxi  •  •  Xi  -  Xf.  (8D) 

4  i  Ai  Ai 

Similarly,  a  sample  value  (x2)y  is  derived  for  the  variable  x,  from  the  integral  equation 

(xp 

J  /2[(x1)y,x2]dx2-/lMOIl]2/,  (8E) 

Ai 

where  A  2  now  has  the  constant  value 

Ai-AJUxOj).  (8F) 

This  procedure  is  continued  until  a  sample  value  has  been  picked  for  each  variable  x.  A  score 
Xy  is  then  calculated  for  the  event  as 

11 


jit 


X,—  C|//l(X|)y . (.X,)j,  .  . .  .(x,)/]. 


(9) 


The  remaining  procedure  is  identical  to  that  already  outlined  in  Subsection  TTA  for 
evaluating  one-dimensional  integrals.  This  procedure  is  outlined  below  for  multqtie  integrals 
and  illustrated  in  Figure  3. 

1.  Pick  a  set  of  sample  values  [(x,  )];  for  the  variables  x,.  These  operations  are  performed 
in  the  following  steps: 

A.  Construct  the  marginal  PDF  /t(x,)  by  using  equation  8A. 

R  Pick  a  sample  value  (x]>,  by  sampling /](X|)  (equation  8R). 

C.  Construct  the  marginal  PDF/2(x2)  by  using  equations  8C  and  8D. 

D.  Pick  a  sample  value  (x2)j  by  sampling  /2(x2)  (equation  8E). 

E  In  a  similar  manner,  pick  a  sample  value  (x,)j  for  each  remaining  x>  by  using  the 
appropriate  marginal  PDF/((x,).  The  reader  should  note  that  the  marginal  PDF  in  each  case  is 
derived  from  a  joint  distribution  of  dimensionality  decremented  by  one  from  the  preceding 
sampling  event. 

2.  Calculate  the  event  score  X/  as  C\  H [(xj)/, ....  (x;  .  ( x,),]. 

3-9.  Calculate  X  and  8X  as  described  earlier  in  this  section.  These  steps  are  identical  to 
those  already  described  in  Subsection  TIA  for  a  one-dimensional  integral. 

Step  9  completes  the  outline  of  a  Monte  Carlo  evaluation  of  multiple  integrals.  The  pro¬ 
cedures  outlined  in  Subsections  ITA  and  TIB  will  be  applied  in  the  following  subsection  to  evalu¬ 
ate  a  multiple  integral  having  the  form  of  the  vulnerability  integral  illustrated  in  equation  1. 


C.  The  Vulnerability  Tiuearal 

The  vulnerability  integral  of  equation  1A  can  often  be  evaluated  by  more  expeditious  cal¬ 
culations  than  those  usually  needed  to  evaluate  multiple  integrals  of  equivalent  dimensionality. 
The  division  of  the  integrand  into  a  source  term,  a  perforation  PDF,  and  a  critical  component 
incapacitation  fit  net  ion,  where  the  source  term  is  a  function  only  of  penetrator  birth  variables, 
will  simplify  the  picking  of  sample  states  for  penetratora.  This  gain  in  calculational  efficiency  is 
even  more  pronounced  for  the  general  vulnerability  equation  where  the  integrals  similar  to 
those  used  in  equation  1  are  sucked  to  represent  the  different  stages  in  a  penetrator  history 
(Reference  5). 

We  will  analyze  the  Monte  Carlo  evaluation  of  the  vulnerability  integral  by  reformulating 
equation  1A.  In  the  new  form, 

X- 

JJ/s('o,w*b0)JJ7/(r,  ,wt,bj  lro,wo,bo,g)0(ri,wi,bi)drj  dwj  rfbi  diodw^dbo, 

W  M  «•  eo  M  M 

"J'J*J*5^ro*wo.*>o)f3l»’o.wo,bo)</rodwodbo,  (10A) 


where 

f>(ro.w0,b0)-J’///(r1,wt,bt  lro,wo,bo,g)OCrt.wi,b|)dr|dw,dbi.  (10B) 

mm  mm  mm 
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Figure  3.  The  Mo  me  Carlo  Estimation  of  a  Multi-Dimensional  Integral 
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If  ZXro.wo.b0)  were  a  known  function  of  the  variables  (ro.wo.bo),  then  X  could  be  evaluated  by 
using  the  Monte  Carlo  method  to  pick  sample  values  (tom  5(ro,wp,bo)  as  outlined  in  Subsec¬ 
tion  TTB.  However,  in  vulnerability  calculations,  ZXro.wo.bo)  is  not  generally  k>  tractable  and  an 
alternate  method  of  solution  must  be  used. 


In  a  Monte  Carlo  solution  still  based  on  the  strong  law  of  large  numbers,  the  multiple 
integrals  in  equation  10A  could  be  evaluated  by  picking  sample  birth  projectiles 
l(ro);,(w0),,(bo)y]  from  the  source  term  5(ro,wo,bo)  by  using  the  procedures  outlined  in  Sub¬ 
section  rm,  and  then  evaluating  each  integral  Z?[(ro)y,(wo)y,(bg)y]  which  is  associated  with  a 
sample  birth  projectile.  The  Monte  Carlo  evaluation  of  each  Z>[(rp)y,(wo)y,(bo)y]  could  also  be 
conducted  by  using  the  procedures  outlined  in  Subsection  TTR  That  is,  a  sufficient  number  K 
of  sample  residual  projectiles  is  picked  from  each  normalized  perforation  PDP 
/^^.wj.b,  t(r0)y,(w0)y,(bo)/«l;A  defined  below,  to  obtain  an  accurate  estimate  Dj  of  the  asso¬ 
ciated  integral  Z?[(r0)y,(w0)y,(bo)y]  where  each  estimate  is  given  by 

Z)y“(Pc)y  (IlA) 

* -1 

The  quantity  X  can  then  be  approximated  as  the  mean  of  a  large  number  of  DJt  that  is 

x=  ^2:  Wj  ZOKnWwjWb,)*].  (11B) 


The  normalized  perforation  PDF  is  given  by 


/^.wi.b,  I (rQ>y,( W0)y,(bo)y ,gj  1 


/fr^wt.b)  l(r0)y,(w0)y,(h0)y,syl 


(11C) 


where  the  normalization  factor  (Pc)y  is  the  probability  that  a  projectile  birthed  as 
[(ro)y,(w0);>(bo).(  )J  w*11  perforate  the  barrier  and  impact  the  critical  component,  that  is 


PcKr0)y.(wo)y.(bo)yI- 

,Wt,b)  1  (rojy.fWojy.Otojy.gyJdr,  rfWj  tfb|.  (llD) 

m  M  aa 


However,  the  approximation  implied  in  equation  1IB  isn’t  necessary  and  the  absolute  con¬ 
vergence  of  X  to  X  can  be  obtained  by  applying  the  Weak  Law  of  Large  Numbers  for  the  case 
where  TchebychefTs  Theorem  is  applicable  (Reference  4).  In  a  Monte  Carlo  estimation  based 
on  this  law,  X  is  estimated  as 

^-7Z^c)z<?[(ri)^W,)y.(b,)yl  (12) 

1 

where  only  one  sample  set  of  values  [(ri)y,(wi)y,(b|  ),1  «  picked  from  each  normalized  perfora¬ 
tion  PDF  /f[ri,W],bi  Kr0)y,(wo)y,(b0)J.  As  noted  in  Reference  4,  the  statement  concerning 
the  convergence  of  the  sample  mean  X  to  the  universe  mean  X  is  weaker  when  only  the  weak 
law  of  large  numbers  is  applicable. 

A  step-by-step  outline  of  a  Monte  Carlo  estimate  of  the  vulnerability  integral  (equation 
IA)  is  outlined  below  and  illustrated  in  Figure  4.  Stacked  vulnerability  integrals  such  as  those 
found  in  vulnerability  methodologies  (Reference  5)  can  be  estimated  by  reiterating  the  follow¬ 
ing  procedure. 

1.  Pick  a  birth  projectile  by  sampling  *S,(r0,w0,b0).  This  sample  projectile  is  identified  as 
[<r0)y,(wo)„(bo)yl.  Since  [r0,wo,bo]  are  each  assumed  to  be  composed  of  several  variables,  the 
procedures  outlined  in  Subsection  ITR  will  have  to  be  used. 
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START 


Figure  ♦.  The  Monte  Carlo  Estimation  of  the  Expected-Value-of-Kill  Vulnerability  Integral 
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2.  Determine  the  values  %t  for  the  variables  associated  with  the  barrier  which  is  impacted 
by  the  projectile. 

3.  Derive  the  perforation  PDF  /[r,tw,,bj  I  (ro)y,(w0)y,(bo)y,|y]  which  predicts  the 
description  (states)  of  the  residual  projectile  at  impact  with  the  critical  component. 

4.  Determine  if  a  residual  penetrator  impacts  the  critical  component  by  comparing  a  ran¬ 
dom  number  RN(0,l)y  with  /,c((ro)y,(wo)y,(bp)y).  Set  the  event  score  to  zero  and  go  to  step  8 
if  a  residual  penetrator  does  not  impact  the  critical  component  This  'Russian  Roulette*  played 
with  the  residual  fragment  differs  from  the  procedures  implied  by  equation  12,  but  the  two 
methods  will  both  converge  toward  the  true  value  of  X. 

5.  Pick  a  residual  peneuator  at  impact  with  the  critical  component  by  sampling  the  nor¬ 
malized  perforation  PDF  /^[rj.wi.bi  l(ro)y,(wo)y,(bo)y,gy].  This  sample  residual  penetrator  is 
identified  as  [(r,)y,(w,)y,(b,  )J. 


6.  Calculate  the  incapacitation  of  the  critical  component  expected  from  the  impact  This 
incapacitation  score  is  identified  as  X;  and  is  given  by 

Xy  -  Q[  (r, ), ,( w,  )y  ,(b,  )j  ].  (13B) 


7.  Accumulate  Xj  and  xj  in  bins  reserved  for  these  operations. 

8.  Conduct  a  total  of  J  similar  events  by  reiterating  Steps  1-7. 

9.  Calculate  an  estimate  of  X  as 

X-^fx,. 

Jj-l 


10.  Calculate  an  estimate  of  8X  as 


*X- 


£x/-/P 

j- j _ 


JU- 1) 


(130 


(13D) 


11.  Assess  8X  to  determine  if  it  is  too  large.  If  necessary,  calculate  more  histories  and 
merge  their  results  with  those  already  calculated  in  order  to  reduce  8X. 

Step  11  completes  the  outline  of  a  Monte  Carlo  evaluation  of  a  multiple  integral  having 
the  form  illustrated  by  the  vulnerability  integral  used  in  equation  1  A.  The  procedures  outlined 
above  will  be  applied  in  the  next  section  to  evaluate  a  definite  integral  constructed  by  using 
simple  analytic  functions. 


ITT.  SAMPLE  PROBLEMS 


a.  A  Monte  Carlo  Solution  Which  Uses  a  General  Multivariate  PDF 

We  have  devised  a  sample  two-dimensional  integral  which  can  be  used  to  illustrate  the 
procedures  outlined  in  both  Sections  ITB  and  ITC.  In  the  selected  definite  integral, 


X  -//  x v  dx  dy  -225/16  - 14.0625, 


(14A) 


I  1 


16 


the  procedures  of  Subsection  ITB  an  illustrated  by  first  rearranging  the  integrand  to  the  form 

(,4B) 

’i  i  * 

Equation  14B  can  then  be  further  changed  to  the  form  used  by  equation  7  A,  that  is 


A —Ci  J*  f H(xj)f(xj)dxdy,  (15A) 

by  taking: 


r 

C,  4’ 

(15B) 

Hbcy)— x2^. 

(150 

and 

(15D) 

when 

(l,l)<(x,y)<(2,2). 

(15E) 

and 

/(xo')-O, 

(15F) 

when 

(xj'XO.l),  (x,y)>(2J). 

(15G) 

The  constant  C \  —9/4  is  introduced  to  normalize /(xj*)  id  unity  over  the  integration  range  of 
the  integral,  that  is, 

/  jf/Ocof )dr ajy -JJ*^  dxdy-l.  (16) 


A  Monte  Carlo  estimate  of  A  can  now  be  calculated  by  picking  sample  values  of  (xo')  from 
/(x,y)  and  then  calculating  the  associated  scores  using  /f(x*y).  The  calculations!  procedure  is 
given  below  and  is  illustrated  in  Hgure  5. 

1.  Pick  a  set  of  values  (x,  j»y)  by  sampling  /(x%y)  over  the  integration  range  of  the  prob¬ 
lem.  This  procedure  is  accomplished  in  the  fallowing  steps 

A.  Construct  the  marginal  PDF/i  (x)  as 

/,<*>- (17a) 
where  f\  (x )  is  taken  to  be  zero  when  x  lies  outside  the  integration  range  ( 1,2). 

R  Pick  a  sample  value  x}  by  solving  the  integral  equation 

X* 

J^<*-/W(0,1)„.  (17B) 

l  3 


C  Construct  a  marginal  PDF/2(x/%y )  fbr  y  as 

„ ,  *  Axj,y) 

/j(W>— 7 - . 

v-2 


(170 
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Figure  3.  The  Mona  Carlo  Estimation  of  the  Example  Problem 
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O.  Pick  a  sample  value  yt  for  y  by  solving  the  integral  equation 


fuxj^dy-RN^Du. 

1 

The  quantity  (Cj)y  can  be  easily  shown  to  be 

so  that  equation  17G  simplifies  to 
f&dy-RNiO,  IV 


(17D) 


(17E) 


C17F) 


(17G) 


1  Calculate  the  score  Ay  for  the  sample  event  as 


(17H) 


1  Calculate  Ay.  Accumulate  Ay  and  k}  in  the  bins  reserved  for  this  operation. 

4  Calculate  a  total  of  J  similar  sample  events  by  reiterating  steps  1-3. 

3.  An  estimate  k  of  k  is  calculated  as 

(171) 

J  j-\ 


6.  An  estimate  of  the  standard  deviation  Ik  of  k  is  calculated  as 


8A- 


j-* 


JU~  1) 


(17J) 


7.  Assess  8A  to  determine  if  it  is  too  large.  If  necessary,  calculate  more  histories  and 
merge  their  results  with  those  already  calculated  in  order  to  reduce  8 A. 

Step  7  completes  an  outline  of  the  calculation  of  an  estimate  of  A  by  using  the  Monte 
Carlo  procedures  given  in  Section  TIB.  A  computer  program  MCTP1  for  performing  this  calcu¬ 
lation  was  written  in  BASIC  and  is  given  in  the  appendix  at  the  end  of  this  report  The  results 
of  calculations  using  MCTP1  are  given  in  TABLE  1  at  the  end  of  this  section 


B.  A  Monte  Carlo  Solution  Which 1 


The  procedures  of  Subsection  TIC  are  illustrated  by  first  rearranging  the  integrand  of  the 
example  integral  (equation  14A)  to  the  form 


i  i 


(ISA) 
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where  the  integration  is  assumed  to  proceed  from  the  inner  integral  to  the  outer  integral.  Equa¬ 
tion  18A  can  be  compared  with  equation  1 A  as 


X- 


f fs(x)/(y  \x)Q{y)dydx -  f}V(x)/*(y  \x)W{x)Q{j)dydx. 

t  t  1*1 

(18B) 

by  taking  ‘  j 

S{x)-x*,  (19A)  j 

0(y)-y.  (19B)  ! 

1 

and 

i 

f(y\x)=xy1.  (19 C)  jj 

The  quantity  W  (x),  identified  here  as  a  weighting  function,  is  introduced  to  compensate  for 
the  normalization  of  the  preceding  source  term  and  PDF  to  unity. 

We  define  a  normalized  source  term  S\x)  as 

St(x)-l^k> 

(19D) 

where 

J3SOf)<fr—  jV(x)dx-/^  <fr-l. 

l  '  i  i  ' 

(19E) 

Applying  the  Theorem  of  Bayes10,  a  normalized  PDF  f*(y  lx)  is  associated  with  xy2  by 

1 

(19F)  J 

or 

/t(ylx)-^-. 

(19G) 

Equation  18A  can  now  be  rewritten  as 

l 

x  -  JJt  sT(x)  Q(y)<fydx , 

i  i  •>  J 

(20a) 

and  rearranged  to 

1 

X=»J*JV(x)/t(,y lx)  Qiy)dydx. 

ii  l3  3  , 

(20B)  - 

I 

The  weighting  function  is  given  by 

FF(x)-y*y  (20C) 


where  —  corresponds  to  the  quantity  Pc  used  earlier  in  outlining  a  solution  of  the  vulnerabil¬ 
ity  integral  (Subsection  IIC). 

i 


Kim,  ■Statistical  Analysis  for  Induction  and  Decision,”  Tha  Drydan  Press,  Inc.,  Hinsdale,  Oinola,  1973, 
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C  END  ") 


Figure  6.  The  Monte  Carlo  Estimation  of  the  Example  Problem  where  the  Integrand  Is  Struc¬ 
tured  in  a  Form  Similar  to  That  of  the  Vulnerability  Integral 


Equation  2 OB  is  in  a  form  which  is  tractable  to  Mona  Carlo  evaluation.  A  general  solution 
is  given  below  and  outlined  in  Figure  6.  The  steps  are: 

1.  Pick  a  sample  value  Xj  for  x  by  sampling  5*  Cx).  The  sample  value  is  derived  from  the 
integral  equation 


(21A) 


7 

The  quantity  y  in  equation  20A  is  saved  and  used  later  as  a  factor  during  the  calculation  of 
the  event  score. 

1  Construct  a  conditional  PDF  f*(y  I xt)  to  be  used  to  pick  a  sample  value  yj  for  y.  The 
conditional  PDF  used  here  (equation  19G)  is  independent  of  Xj,  but  this  independence  may  not 
exist  for  cases  where  functions  other  than  xy1  are  used  in  the  construction. 

3.  Pick  a  sample  value  yt  by  sampling  f*(y  I  Xj).  This  sample  value  is  derived  from  the 
integral  equation 


j/V*;  )dy  -  jN?i  dy -[/W((U>H/. 

1  l  ' 


lx 

The  quantity  ■—  in  equation  20A  is  saved  and  used  later  as  a  factor  during  the  calculation  of 
the  event  score. 


4.  Calculate  the  score  for  the  event  as 

5.  Calculate  x/.  Accumulate  X,  and  X,J  in  the  bins  reserved  for  this  operation. 


(21C) 


6.  An  estimate  X  of  X  is  calculated  as 
_  i  j 
tZV 

J  j- t 


(21D) 


7.  An  estimate  8X  of  the  standard  deviation  of  X  is  calculated  as 

H^iT 


8.  Assess  5X  to  dearmlne  if  it  is  too  large.  If  necessary,  conduct  more  histories  by 
reiterating  Steps  1-5  and  merge  these  results  with  those  obtained  earlier. 

Step  S  completes  the  outline  of  the  Monte  Carlo  evaluation  of  the  example  problem  by 
using  the  procedures  which  would  be  used  to  evaluate  a  vulnerability  integral  of  the  form 
shown  in  equation  1.  The  results  of  test  calculations  of  this  example  problem  are  given  below 
in  Subsection  TTTC. 
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C  Example  Problem  Results 


A  computer  program  was  written  in  BASIC  (APPENDIX)  to  solve  the  sample  integral  by 
using  both  Monte  Carlo  procedures  described  in  the  preceding  sections.  The  results  of  test  cal¬ 
culations  where  5000  histories  were  conducted  Tor  each  estimate  are  given  below  in  Table  I. 
The  reader  should  note  that  standard  deviations  of  magnitude  ten  times  that  of  the  sample  cal¬ 
culations  could  be  obtained  by  using  approximately  fifty  histories.  Confidence  levels  of  this 
order  of  magnitude  would  usually  be  acceptable  in  vulnerability  calculations 


IV.  CONCLUSION 


We  have  outlined  the  use  of  the  Monte  Carlo  method  for  solving  expected  value  integrals 
of  the  type  encountered  in  ballistic  vulnerability  calculations.  As  an  aid  to  the  comprehension  of 
these  methods,  the  outlined  procedures  (Subsection  TTB  and  ITC)  were  applied  to  a  sample 
problem  for  which  a  closed  form  solution  existed  and  could  be  readily  calculated.  The 
insignificant  differences  between  the  computer -generated  Monte  Carlo  solutions  and  the  exact 
solution  (Table  1)  serve  to  indicate  the  viability  of  the  Monte  Carlo  techniques. 

The  reasons  for  using  the  Monte  Carlo  method  in  preference  to  deterministic  methods  to 
solve  vulnerability  problems  can  be  summarized  as. 


1.  Regardless  of  the  dimensionality  af  die  integration,  a  Monte  Carlo  estimate  of  multiple 
integrals  converges  toward  the  true  value  as  -J=-  where  N  is  the  number  of  trials  used  in  calcu¬ 


lating  the  estimate.  Therefore,  the  Monte  Carlo  method  may  be  the  most  efficient  procedure 
for  solving  expected -value  problems  when  the  dimensionality  of  the  associated  multiple 
integrals  is  large  and  confidence  levels  for  the  estimate  on  the  order  of  five  percent  are  accept¬ 
able. 


2.  The  probability  of  introducing  systematic  errors  such  as  those  which  might  be  intro¬ 
duced  in  deterministic  calculations  using  a  computerized  phantom  of  the  target  vehicle  is 
greatly  reduced. 

3.  Calculational  simplicity  can  often  be  obtained  by  the  division  of  the  integrand  into  one 
part  used  in  calculating  scores  and  another  pan  used  to  derive  the  sampling  PDFs.  In  particu¬ 
lar,  the  PDFs  used  in  describing  the  ballistic  phenomena  are  often  directly  derivable  in  terms 
of  the  "used"  random  variables  B  and  G. 

4.  The  calculations  can  bs  organized  so  that  the  sample  events  correspond  in  a  one-to- 
one  manner  with  actual  ballistic  events.  This  organization  would  aid  vulnerability  analysts  in 
conducting  calculations  and  interpreting  results. 

5.  In  many  vulnerability  problems,  sample  events  can  be  obtained  more  easily  than  the 
associated  PDF.  Such  problems  must  generally  be  solved  by  using  the  Monte  Carlo  method. 
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METHOD  OF 
SOLUTION 

RUN 

NO. 

ESTIMATE 

STANDARD 

DEVIATION 

VULNERABILITY 

INTEGRAL 

1 

14.058 

0.048 

2 

14.158 

0.048 

3 

14.043 

0.048 

MUITWARIATE 

SAMPLING 

1 

14.153 

0.102 

2 

14.082 

0.101 

3 

14.044 

0.100 

CLOSED  FORM 

14.0625 

Table  1.  The  Mome  Carlo  Estimates  of  the  Example-Problem  Integral 
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appendix,  the  computer  program  used  to  solve  the  example 

PROBLEMS 


1  REM  MC  TEST  PROG  1  MCTP1 
10  INPUT  "PI,  N"  ■  PI,  N 
iO  FOR  I=1T0  N 
10  GOSUB  70 
40  X*SQR<1+2*R> 

50  GOSIJB  '  70 
6 0  ysbOP.:  1+3 ♦  p 
70  2=X*V+X*V 
SO  S1=S1+Z 
SO  S2=S2+Z+Z 

100  SELECT  PRINT  005:  PRINT  I,  Z 

110  NEXT  I 

120  S1=S1XN 

ICO  S2=S2-N*S1*S1 

140  S2=S2/ N*  <  N-l 

150  S2=S0P<S2) 

ISO  SELECT  PRINT  215 
l'O  PRINT  9*51/4, 3+S2/4 
ISO  END 

4000  OEFFN  •"  70 

4001  REM  RN  SP 
40i05  R1=25*R1 
4 Ovio  GOSUB  '71 
4015  r-  i=5*Rl 

4 0 2 0  GOSIJB  ’71 
4025  F=P1.’S710S864 

4o:o  return 

4050  OEFFN  •'  71 

4051  REM  RN  SP2 
4055  R2=R1/S71088S4 
40S0  R3  =S7108SS4* 1  NT R2  5 
4 OS 5  R1-P1-R2 

4070  RETURN 


1  REM  MC  TEST  PROG  2  MCTP2 

10  INPUT  "Rl,  N",  PI,  N 

20  FOR  I-1T0  N 

30  ijQSlJE'  70 

40  X=*:  1+7+R  '  ~0.  333333 

50  GOSUB  '70 

60  Vs <  1+7  +R  j  "0.  233222 

70  Z=X+V 

30  Sl-Sl+Z 

30  S2=S2+Z+Z 

1O0  SELECT  PRINT  005:  PRINT  I 

HO  NEXT  I 

120  S1=S1/N 

120  S2=S2-N*S1*S1 

140  S2=S2/ <  N* <  N-l > > 

150  S2=S0R  •.  S2  > 

ISO  SELECT  PRINT  215 
170  PRINT  49+S1/9,  49+S2/9 
130  END 

4000  OEFFN '70 

4001  REM  RN  SP 
40O5  R1=25*R1 
4010  GOSUB  -'71 
4015  R1=5*R1 
4020  GOSUB  71 
4025  R=R1.‘671088S4 
4030  RETURN 

4050  OEFFN "71 

4051  REM  RN  SR2 
4055  R2=R1 X S71083S4 
40SO  R3=S7 108864* I NT < R2  J 
40S5  P1=R1-R3 

4070  RETURN 
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